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ABSTRACT 

We present a new algorithm for identifying dark matter halos, substructure, and tidal features. The approach 
is based on adaptive hierarchical refinement of friends -of-friends groups in six phase-space dimensions and one 
time dimension, which allows for robust (grid-independent, shape-independent, and noise-resilient) tracking of 
substructure; as such, it is named Rockstar (Robust Overdensity Calculation using K-Space Topologically 
Adaptive Refinement). Our method is massively parallel (up to 10^ CPUs) and runs on the very largest simu- 
lations (>10'" particles) with high effi ciency (10 CPU h ours and 60 gigabytes of memory required per billion 
particles analyzed). A previous paper (|Knebe et al.|201 1) has shown Rockstar to have class-leading recovery of 
halo properties; we expand on these comparisons with more tests and higher-resolution simulations. We show a 
significant improvement in substructure recovery as compared to other halo finders and discuss the theoretical 
and practical limits of simulations in this regard. Finally, we present results which demonstrate conclusively 
that dark matter halo cores are not at rest relative to the halo bulk or satellite average velocities and have co- 
herent velocity offsets across a wide range of halo masses and redshifts. For massive clusters, these offsets 
can be up to 400 km s~' at z = and even higher at high redshifts. Our implementation is pubUcly available at 
http: //code. google. com/p/ rockstar 

Subject headings: dark matter — galaxies: abundances — galaxies: evolution — methods: N-body simulations 



1. INTRODUCTION 

In the paradigm of Lambda Cold Dark Matter (ACDM), 
the majority of the matter density of the universe does not 
couple to electromagnetic fields, leaving it detectable only 
through its gravitational and possibly weak force interactions. 
Nonetheless, the effects of dark matter on the visible universe 
are spectacular, as the steep gravitational potentials of bound 
dark matter halos channel baryons together, forming the birth- 
places of visible galaxies. In this model, the locations, sizes, 
and merging histories for galaxies are thus intricately con- 
nected to the growth of bound dark matter structures. 

Testing this model in detail requires extensive computer 
simulations, as the complicated nonlinear evolution of struc- 
ture growth cannot be fully evaluated by hand. However, 
the simulations, which in general track the evolution of a 
set of dark matter particles, must be postprocessed to de- 
termine the locations and properties of bound dark matter 
structures also known as "halos" — namely, the locations and 
properties which influence the formation of visible galaxies. 
This postprocessing necessarily involves both ambiguity and 
imprecision — ambiguity in the definitions (e.g., the center of 
a bound halo) and imprecision in determining halo proper- 
ties due to limited information (e.g., for halos consisting of 
only a few dark matter particles, or for determining particle 
membership in two overlapping halos). Currently, as a con- 
sequence, essential statistical measures (e.g., the number den- 
sity of halos with a given mass) are known to at best 5-10% 
even for a specific cosmology ( [Tinker et al.||200 8). More- 
over, halo properties have rarely been tested for consistency 
between timesteps, which can be a serious problem for robust 
modeling of their historie s required for any theory of galaxy 
forma tion (see, however, [Behroozi et al.|[2011[ [Tweed et al.| 
|2009] l. 

This level of uncertainty is unacceptable for future surveys, 
including data expected to come from DES, Pan-STARRS, 



eROSITA, Herschel, Planck, JWST, and LSST, in order to 
fully realize the constraining power of these observations. In- 
deed, derived quantities such as the halo mass function and 
autocorrelation functions must be understood at the one per- 
cent level in order for theoretical uncertainties to be at the 
same level as statisti cal unce rtainties in constraining, e.g., 
dark energy ( [Wu et al 2010; Cunha & Evrard 2010 for exam- 
ple). While some of the current uncertainties may be pinned 
on dark matter simulations ( e.g., the effects of baryons on halo 
profiles; Stanek et al.|2009 1, different halo finders running on 
the same simulation demonstrate that a large fraction is s till 
due to the process of halo finding itself ( jKnebe et al.|201 1} . 

As previously mentioned, some of these uncertamties are 
due to limited use of information. Considerable progress has 
been made since the first generation of position-space halo 



finders ( [Davis et al.||1985|]Lacey & Cole|[1994) , which used 
only particle locations to determine bound halos. Currently, 
the most advanced al gorithms are adaptive phase-s pace find- 
ers (e.g., Maciejewsk i et al.||2009| [Ascasibar et al. 20101), 
which make use of the full six dimensions of particle positions 
and velocities. At the same time, an often-overlooked (with 



the possible exception of Tweed et al. 2009j ) aspect of halo 
finding is the extra information stored in the temporal evolu- 
tion of bound particles. In this paper, we detail an advanced 
adaptive phase-space and temporal finder which is designed 
to maximize consistency of halo properties across timesteps, 
rather than just within a si ngle simulation tim estep. Together 
with a companion paper ( Behroozi et al.poTi) , which details 
the process of comparing and validating halo catalogs across 
timesteps to create gravitationally consistent merger trees, our 
combined approach is the first to use particle information in 
seven dimensions to determine halo catalogs, allowing un- 
precedented accuracy and precision in determining halo prop- 
erties. 

Furthermore, in contrast to previous grid-based adaptive 
phase-space approaches, ours is the first grid-independent and 
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orientation-independent approach; it is also the first publicly- 
available adaptive phase-space code designed to run on mas- 
sively parallel systems on the very large simulations (> 10'° 
particles) which are necessary to constrain structure formation 
across the range of scales probed by current and future galaxy 
observations. Finally, we remark that our halo finder is the 
first which is successfully able to probe substructure masses 



down to the very centers of host halos (see also Knebe et al. 
201 1\ , which is essential for a full modeling of galaxy statis- 



tics and will enable future studies of the expected breakdown 
between halo positions and galaxy positions due to the effects 
of baryon interactions at the very centers of galaxies. 

Throughout, we have paid careful attention not only to the 
basic task of assigning particles to halos, but also to the pro- 
cess of estimating useful properties from them to compare 
with observations. While in many cases galaxy surveys are 
not able to probe halo properties to the same precision as halo 
finders in simulations, one significant counterexample exists. 
It is a com mon practice especially for FOF-based halo find- 
ers (e.g., Davis et al.'1985','Springel et al.'200T]|Habib et al. 



12009 ; Rasera et al. 2010 ; Gardner et al. 2007) to calculate halo 
velocities by averaging all halo particle velocities together to 
find a bulk velocity. Examination of the difference between 
velocities in the inner regions of halos and the bulk averaged 
velocity suggests that the bulk average velocity may be off- 
set by several hundred km s"' from the velocity at the radius 
where galaxies are expected to reside; differences at this scale 
are easily detectable in cluster redshift surveys. As this dif- 
ference has an important impact on the usefulness of derived 
halo properties, we additionally perform an investigation of 
the core-bulk velocity difference in halos across a wide range 
of redshifts and masses. 

We begin this paper with a survey of previous work in halo 
finding as well as previous limitations in ^ We discuss our 
improved methodology in ^and conduct detailed tests of our 
approach in ^ We present an analysis of the theoretical and 
practical limitations of simulations in terms of tracking sub- 
structure in ^ Finally, our results concerning the velocity 
offsets of cluster cores are presented in ^ we summarize 
our conclusions in ^ Multiple simulations including slight 
variations of cosmological parameters are considered in this 
paper; all simulations model a flat ACDM universe and we 
always take the Hubble constant Hq to be 70 km s"' Mpc~'; 
equivalently, h = Q.7. 

2. PREVIOUS HALO FINDERS 

2.1. Summary of Previous Approaches 

Previously published approaches to halo finding may be 
classified, with few exceptions, into two large groups. Spher- 
ical overdensity finders, such as AHF, ASOHF, BDM, 
SO, pSO, VOBOZ, and SKI D (|Knollma nn & Knebe||2009t 
Planelles & Quilis 2010, Kly pin et al. p'999; Lacey & Cole" 
B94, Jenkins et al. 2001, Evrard et al.l2002; Sutter & Ricker, 



2010} [Neyrinck et al.|20()5| [Stadel 2001) proceed by identify- 
ing density peaks in the particle distribution and then adding 
particles in spheres of increasing size around the peaks until 
the enclosed mass falls below a predetermined density thresh- 
old (a top-down approach). Friends-of-friends and HOP- 
based halo finders, such as EO F, SUBFIND LANL, pFOF, 
Ntropy-fofsv, and AdaptaHOP (|Davis et aL]|1985 ', 'Springell 
eTarillOOT] |Habib et al.|[2009r|kasera et al.| |^010; Gard- 



designed, search for substructure inside these particle collec- 
tions (a bottom-up approach). Phase-space finders typically 
extend these two approaches to include particle velocity in- 
formation, either by c alculating a phase-space density, such 
as HSF and HOT6D ( [Maciejewski et al ."2009; Ascasibar etj 
,al. ,2010) or by using a phase -space linking length, as does 
WVOV piemand et al.|2006l ). 

There are two notable exceptions to these algorithms. 



ner et al.|2007| [Tweed et all 
which fall above a certain at 



[2009 1, collect particles together 



ensity threshold and then, if so 



name ly the ORI GAMI halo finder (discussed in Knebe et al. 
[20n] l and HBT ( |Han et al. l [20n] l. ORIGAMI operates by 
examining phase-space shell crossings for the current parti- 
cle distribution as compared to the initial particle distribution; 
shells which have crossed along three dimensions are con- 
sidered to be halos (as opposed to shells which have crossed 
along one or two dimensions, which would be considered as 
sheets and filaments, respectively). HBT uses a friends-of- 
friends approach to find distinct halos and uses particle lists 
from distinct halos at previous timesteps to test for the pres- 
ence of subhalos. These algorithms both rely heavily on tem- 
poral information in their approach to halo finding. 

2.2. Limitations of Previous Algorithms 

In order to develop an improved halo finder, it is impor- 
tant to understand some of the shortcomings of previous ap- 
proaches. For the vast majority of halos, even the most basic 
of algorithms (EOF and SO) do an acce ptable job of deter - 
mining halo properties to 10% accuracy ( [Knebe et al.|201 1} . 
However, recent interest in the detailed properties and his- 
tories of halos — e.g., precision mass functions and merger 
trees — as well as interest in the shape of tidal structures re- 
quires improvements to older approaches; this has resulted in 
a proli feration of new codes in the past decade (summarized 
in Kne be et al.|201 l| l. 

The most significant improvements to halo finding have 
come from using the information from six-dimensional (po- 
sition and velocity) phase space. Two traditional weak points 
for 3D (position-space) halo finders have been major merg- 
ers and subhalos close to the centers of their host halos. In 
both cases, the density contrast alone is not enough to dis- 
tinguish which particles belong to which halo: when two ha- 
los are close enough, the assignment of particles to halos be- 
comes essentially random in the overlap region. However, as 
long as the two halos have relative motion, six-dimensional 
halo finders can use particles' velocity space information to 
very effectively determine particle-halo membership. This, 
coupled with the ability of 6D halo finders to find tidal rem- 
nants (which are condensed in phase space but not in position 
space), means that phase-space capabilities are required for 
the most accurate and interesting studies of dark matter halos. 

At the same time, phase space presents a unique challenge. 
While position space and velocity space have well-defined 
distance metrics, there is not a single, unique way to com- 
bine position and velocity distance into a single metric. For 
a useful definition of phase-space distance, one needs to be 
able to decide, e.g., whether an object just passing by the ori- 
gin at 1 km s~^ is closer or farther than an object at rest 1 
kpc away. One approach, used by 6DFOF, is to choose in ad- 
vance a static conversion between velocity and position space. 
While simple, this approach seems both somewhat arbitrary 
(why choose a particular value for the conversion?) and self- 
limiting: if too short a linking-length is chosen, the full ex- 
tent of substructures cannot be found; if too large a linking- 
length is chosen, then otherwise distinct substructures will be 
merged together. 
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A demonstrably superi or approach, at lea st in terms of re- 
covering halo properties ( jKnebe et al.|201l} , is to adoptively 
define a phase-space metric. Both HSF and HOT6D subdi- 
vide the simulation space into six-dimensional hyperboxes 
containing (at the maximum refinement level) as little as a 
single particle each. For a given particle, the enclosing hy- 
perbox gives a local estimate of the phase-space metric, based 
on the relative sizes of the hyperbox's dimensions in position 
and velocity space. The usefulness of this estimate depends 
heavily on the method for partitioning space into hyperboxes; 
HSF uses, for example, an entropy-based approach to deter- 
mine whether more information is contained in the particle 
locations for position space or velocity space. 

These algorithms all give excellent results for identifying 
halo centers at a single timestep. However, consistent halo 
catalogs across timesteps are often compromised by a funda- 
mental ambiguity in the definition of a host halo. For major 
mergers, it is often unclear which halo center represents the 
larger "host" or central halo, and which represents the sub- 
halo. Phase-space halo finding helps when the two halo cen- 
ters are relatively far apart (i.e., weakly interacting), because 
there exists a strong correlation between the velocities of par- 
ticles in the halo outskirts and halo centers. However, when 
the centers come close enough to interact strongly, this cor- 
relation is weakened, and it becomes much more difficult to 
accurately assign particles to the halos. As a result, it is much 
more difficult to determine which of the halo centers should 
be considered the host halo. Since the definition of halo mass 
often includes the mass of subhalos, this problem can result 
in large mass fluctuations across timesteps for merging halos. 
A number of solutions to this problem have been pro- 



posed and ex amined with the Ada ptaHOP halo finder (Tweed 
1^ al. 2009). [Tweed et af] ( [2009] l found that a temporal ap- 



proach (examining the host vs. subhalo assignment at earlier 
timesteps) was most successful at fixing this problem. Other 
methods, such as choosing the densest halo center to be the 
host halo, have inherent instabilities because of the spread in 
halo concentrations at a g iven halo mass. At the same time, 
while Tweed et al. ( |2009| l successfully resolves this problem, 
it nonetheless only finds halos in position space, and thus has 
the same weaknesses in identifying subhalo and major merger 
properties. 



The HBT algorithm ( |Han et aL 2011 1 uses the ingenious 
approach of tracing subhalos by using the particles found in 
previously distinct friends-of-friends halos, which could po- 
tentially also solve many of these problems. Yet, it is not 
perfect because it also includes assumptions about accretion 
onto subhalos (e.g., that subhalos cannot accrete background 
matter from the host) which are untrue especially for large 
linking-lengths (as halos will be identified as satellites far 
from the actual virial radius of the host) and major mergers 
(where satellite and host are more ambiguous; they can in any 
case easily trade particles with each other). These vastly sim- 
plify the code at some expense to the completeness and ac- 
curacy of the mass function. More seriously, the design of 
the algorithm requires temporal information to find subhalos; 
in cases where simulation outputs are spaced very far apart 
or when only one timestep is available, it cannot effectively 
function to find substructure. These issues are fixable: a fu- 
ture version of the algorithm could easily combine advanced 
single-timestep substructure finding with checks against pre- 
vious timesteps' particle lists; however, future halo finders are 
not the subject of this paper. 



3. AN IMPROVED APPROACH: ROCKSTAR 

Our primary motivation in developing a halo finder was to 
improve the accuracy of halo merger trees that are required 
for an understanding of galaxy evolution. The design of our 
halo finder was thus motivated by a requirement for consis- 
tent accuracy across multiple timesteps. This interest led to 
the development of a unique, adaptive phase-space temporal 
algorithm which is provably independent of halo orientation 
and velocity relative to the simulation axes, and which also 
attempts to be highly preserving of particle-halo and halo- 
subhalo membership across timesteps. In addition, we paid 
special attention to the algorithm's efficiency and paralleliz- 
ability, to allow it to run on the largest current datasets and 
so that it could easily scale to the next generation of simu- 
lations. Thus far, we have run the halo finder (and in many 
cases the partner merger tre e code) on the Bolshoi (^ 10'° 
particles, Klypin et al.|201 1 1 and LasDamas simulations (200 
boxes of 1 - 4 X lO'^ particles each, McBride et al. in prepara- 
tion), on several 2048"^ simulations run to create Dark Energy 
Survey (DES) simulated sky catalogs, on ^ one hundred high 
resolution halos simulated as part of the RHAPSODY project 
(Wu et al, in preparation), and on halos A-1 through A-5 of 
the Aquari us Simulati on (up to 1.4 x 10'^ particles in a single 
halo; jSpringel et al.|2008) . 

As a first step, our algorithm performs a rapid variant of 
the 3D friends-of-friends (FOF) method to find overdense re- 
gions which are then distributed among processors for anal- 
ysis (§3.1 1. Then, it builds a hierarchy of FOF subgroups 



in phase space by progressively and adaptively reducing the 
6D linking length, so that a tunable fraction of particles are 
captured at e ach subgroup as compared to the immediate par- 
ent group (5 3.2 1. Next, it converts this hierarchy of FOF 
groups into a list of particle memberships for halos (5 3.3 1. 
It then computes the host halo/subhalo relationships among 
halos, using information from the previous timestep if avail- 
able ( ^3.4| i. Finally, it removes unbound particles from halos 
and calculates halo properties, before automatically generat- 
ing particle-based merger trees (§3.5 i. A visual summary of 
these steps is shown in Fig.[T] 

3.1. Efficient, Parallel FOF Group Calculation 

The 3D friends-o f-friends (FOF) a lgorithm has existed 
since at least 1985 ( Davis et al.||1985 i. In principle, imple- 
mentation of the algorithm is straightforward: one attaches 
two particles to the same group if they are within a prespec- 
ified linking length. Typically, this linking length is chosen 
in terms of a fraction b of the mean interparticle distance (in 
our code, as for others, the cube root of the mean particle vol- 
ume); common val ues for generating halo catalogs range from 
b = 0.l5tob = 0.2 (More et al.|201 l| l. 

In practice, this means that one must determine neighbors 
for every particle within a sphere of radius equal to the link- 
ing length. Even with an efficient tree code (we use a cus- 
tom binary space partitioning tree), this represents a great deal 
of wasted computation, especially in dense cluster cores. In 
such cases, particles might have tens of thousands of neigh- 
bors within a linking length, all of which will eventually end 
up in the same FOF group. 

As the 3D FOF groups are used in our method only to di- 
vide up the simulation volume into manageable work units, 
we instead make use of a modified algorithm which is an or- 
der of magnitude faster. As is usual, neighboring particles are 
assigned to be in the same group if their distance is within the 
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1. The simulation volume is divided 
into 3D Friends-of-Friends groups 
for easy parallelization. 

2. For each group, particle positions 
and velocities are divided (normal- 
ized) by the group position and ve- 
locity dispersions, giving a natural 
phase-space metric. 

3. A phase-space linking length is 
adaptively chosen such that 70% of 
the group's particles are linked to- 
gether in subgroups. 

4. The process repeats for each 
subgroup: renormalization, a new 
linking-length, and a new level of 
substructure calculated. 

5. Once all levels of substructure are 
found, seed halos are placed at the 
lowest substructure levels and par- 
ticles are assigned hierarchically to 
the closest seed halo in phase space. 

6. Once particles have been as- 
signed to halos, unbound particles 
are removed and halo properties 
(positions, velocities, etc.) are 
calculated. 

Fig. 1 . — A visual summary of the particle-halo assignment algorithm. 



linking length. However, if a particle has more than a certain 
number of neighbors within the linking length (16, in our ver- 
sion), then the neighbor-finding process for those neighboring 
particles is skipped. Instead, neighbors for the original parti- 
cle out to twice the original linking length are calculated. If 
any of those particles belong to another FOF group, that cor- 
responding group is joined with that of the original particle. 
Thus, although the neighbor-finding process has been skipped 
for the nearest particles, groups which would have been linked 
through those intermediate neighbors are still joined together. 

This process is provably more inclusive than the standard 
FOF algorithm — fine for our desired purpose — as well as 
much faster: neighbors must be calculated over a larger dis- 
tance, but many fewer of those calculations must be per- 
formed. Indeed, rather than slowing down as the linking 
length is increased, our variation of FOF becomes faster Be- 
cause of this, we are free to choose an exceptionally large 
value for the linking length. A value of b = 0.2 is too small, as 
it does not include all particles out to the virial radius (More 




densities can be determined for even the most ellipsoidal ha- 
los. 

As concerns parallelization, the most important work oc- 
curs at this stage. Separate reader tasks load particles from 
snapshot files. Depending on the number of available CPUs 
for analysis, a master process divides the simulation region 
into rectangular boundaries, and it directs the reader tasks to 
send particles within those boundaries to the appropriate anal- 
ysis tasks, along with particles in surrounding "ghost" zones, 
so that properties for halos near analysis boundaries are ap- 
propriately recovered. The recommended and default zone 
width is 3.0 Mpc ; assu ming a halo defintion based on Avi, 
(Bryan & Norman 1998| l, even the largest halos have virial 
radii less than 2.5 Mpc 

Each analysis task first calculates 3D FOFs in its assigned 
analysis region; the FOF groups are then distributed for fur- 
ther phase-space analysis according to individual processor 
load. The load-balancing procedure is described in further 
detail in Appendix [A] In its current version, no group "stitch- 
ing" is performed — groups must lie entirely within single 
analysis regions (including ghost zones), making this version 
of the code less appropriate for parallel use with simulations 
less than 15 Mpc /;"' across. In such cases (e.g., running on a 
single high resolution halo), the code can be run in serial. In 
the current implementation, there is no support for multiple 
particle masses, although support could easily be added by 
varying the linking length depending on particle mass. Pro- 
vided enough interest, these features may be added in a future 
version of Rockstar. 

3.2. The Phase-Space FOF Hierarchy 

For each 3D FOF group which is created in the previous 
step, the algorithm proceeds by building a hierarchy of FOF 
subgroups in phase space. Deeper levels of subgroups have a 
tighter linking-length criterion in phase space, which means 
that deeper levels correspond to increasingly tighter isoden- 
sity contours around peaks in the phase-space density distri- 
bution. This enables an easy way to distinguish separate sub- 
structures — above some threshold phase-space density, their 
particle distributions must be distinct in phase space; other- 
wise, it would be difficult to justify the separation into differ- 
ent structures!]] 

Beginning with a base FOF group, Rockstar adaptively 
chooses a phase-space linking length based on the standard 
deviations of the particle distribution in position and velocity 
space. That is, for two particles pi and p2 in the base group, 
the phase-space distance metric is defined as: 



d{pup2) = 



1/2 



(1) 



et 



al.|201 1 1; after evaluating different choices for b (see ^4.5 
we chose b = 0.28, which guarantees that virial spherical over 



where and (t,. are the particle position and velocity disper- 
sions for the given FOF group; this is identical to the metric of 
Gottlober ( 1998 ). For each particle, the distance to the near- 
est neighbor is computed; the phase-space linking length is 
then chosen such that a constant fraction / of the particles are 
linked together with at least one other particle. In large groups 
(> 10,000 particles), where computing the nearest neighbor for 
all particles can be very costly, the nearest neighbors are only 
calculated for a random 10,000-particle subset of the group, 

' This would not be true if, e.g., halos had Plummer profiles or otherwise 
flat density profiles in their centers. 
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as this is sufficient to determine the Hnking length to reason- 
able precision. 

The proper choice of / is constrained by two considera- 
tions. If one chooses too large a value (/ > 0.9), the algo- 
rithm will take much longer, and it can also find spurious (not 
statistically significant) subgroups. If one chooses too low of 
a value (/ < 0.5), the algorithm may not find smaller sub- 
structures. As such, we use an intermediate value (/ = 0.7); 
with the recommended minimum threshold for halo p articles 
(20). In our tests of the mock NFW halos described in |Knebe| 
|et al.| (j201 1 ), this results in a false positive rate of 10% for 
20-particle groups compared to a cosmological subhalo dis- 
tribution, which declines exponentially for larger group sizes. 
These false positives are easily rem oved by the sig nificance 
and boundedness tests described in 33.3| and §3.5.3 

Once subgroups have been found in the base FOP group, 
this process is repeated. For each subgroup, the phase-space 
metric is recalculated, and a new linking-length is selected 
such that a fraction / of the subgroup's particles are linked 
together into sub-subgroups. Group finding proceeds hierar- 
chically in phase space until a predetermined minimum num- 
ber of particles remain at the deepest level of the hierarchy. 
Here we set this minimum number to 10 particles, although 
halo properties are not robust approaching this minimum. 

3.3. Converting FOF Subgroups to Halos 

For each of the subgroups at the deepest level of the FOF hi- 
erarchy (corresponding to the local phase-space density max- 
ima), a seed halo is generated. The algorithm then recursively 
analyzes higher levels of the hierarchy to assign particles to 
these seed halos until all particles in the original FOF group 
have been assigned. To prevent cases where noise gives rise to 
duplicated seed halos, we automatically calculate the Poisson 
uncertainty in seed halo positions and velocities, and merge 
the two seed halos if their positions and velocities are within 
10(7 of the uncertainties. Specifically, the uncertainties are 
calculated as /iv = (Jx/^/n and /i,, = a^.j ^Jn, where and ct,, 
are the position and velocity dispersions, and n is the number 
of particles in the seed halo. The two halos are merged if 

y^(^i-^2)W + (vi-V2)V;' < 10\/2. (2) 

For a parent group which contains only a single seed halo, 
all the particles in the group are assigned to the single seed. 
For a parent group which contains multiple seed halos, how- 
ever, particles in the group are assigned to the closest seed 
halo in phase space. In this case, the phase-space metric is 
set by the seed halo properties, so that the distance between a 
halo h and a particle p is defined as: 

diKp^=[\^.\^\'" (3) 

where Uy is the seed halo's current velocity dispersion, M/, is 
its current mass, and "vir" spe cifies the virial ov erdensity, us- 
ing the definition of /3i,,> from Bryan & Norman| ( [T998) , which 
corresponds to 360 times the background density at z = 0, 
however, other choices of this density can easily be applied. 

Using the halo radius ^vir as the position-space distance nor- 
malization may seem unusual at first, but the natural alter- 
native (using (Jx) gives unstable and nonintuitive results. At 



fixed phase-space density, subhalos and tidal streams (which 
have lower velocity dispersions than the host halo) will have 
larger position-space dispersions than the host halo. Thus, if 
Ox were used, particles in the outskirts of a halo could be eas- 
ily mis-assigned to a subhalo instead of the host halo. Using 
rvii, on the other hand, prevents this problem by ensuring that 
particles assigned to subhalos cannot be too far from the main 
density peak even if they are close in velocity space|^ Intu- 
itively, the largest effect of using r^„ is that velocity-space 
information becomes the dominant method of distinguishing 
particle membership when two halos are within each other's 
virial radii 

This process of particle assignment assures that substruc- 
ture masses are calculated correctly independently of the 
choice of /, the fraction of particles present in each subgroup 
relative to its parent group. In addition, for a subhalo close to 
the center of its host halo, it assures that host particles are not 
mis-assigned to the subhalo — the central particles of the host 
will naturally be closer in phase space to the true host center 
than they are to the subhalo's center 

3.4. Calculating Substructure Membership 

In addition to calculating particle-halo membership, it is 
also necessary to determine which halos are substructures of 
other halos. The most common definition of substructure is 
a bound halo contained within another, larger halo. Yet, as 
halo masses are commonly defined to include substructure, 
the question of which of two halos is the largest (and thus, 
which should be called a satellite of the other) can change de- 
pending on which substructures have been assigned to them. 
This is one of the largest sources of ambiguity between spher- 
ical overdensity halo finders, even those which limit them- 
selves to distinct halos. 

We break the self-circularity by assigning satellite member- 
ship based on phase-space distances before calculating halo 
masses. Treating each halo center like a particle, we use the 
same metric as Eq. [3] and calculate the distance to all other 
halos with larger numbers of assigned particles. The halo in 
question is then assigned to be a subhalo of the closest such 
halo, if one exists. If the halo catalog at an earlier timestep 
is available, this assignment is modified to include temporal 
information. Halo cores at the current timestep are associ- 
ated with halos at the previous timestep by finding the halo at 
the previous timestep with the largest contribution to the cur- 
rent halo core's particle membership. Then, host-subhalo re- 
lationships are checked against the previous timestep; if nec- 
essary, the choice of which halo is the host may be switched 
so as to preserve the host-subhalo relationship of the previous 
timestep. 

As explained above, these host-subhalo relationships are 
only used internally for calculating masses: particles assigned 
to the host are not counted within the mass of the subhalo, but 
particles within the subhalo are counted as part of the mass 
of the host halo. This choice assures that the mass of a halo 
won't suddenly change as it crosses the virial radius of a larger 
halo, and it provides more stable mass definitions in major 
mergers. Once halo masses have been calculated, the subhalo 
membership is recalulated according to the standard definition 

^ For determination of tidal streams, this "problem" becomes a "feature," 
and use of CTv may be preferable to r^\^. 

' An alternate radius (e.g., r-itxib or ''soor) could be used instead, but it 
would have an effect only on a small fraction of particles in a small fraction 
of halos (major mergers). 
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(subhalos are within ta of more massive host halos) when the 
merger trees are constructed. 

For clarity, it should be noted that every density peak within 
the original FOF analysis group will correspond to either a 
host halo or a subhalo in the final catalog. It has been observed 
that FOF groups will "bridge" or "premerge" long before their 



corresponding SO halo counterparts (e.g., Klypin et al. 20lT| l. 
However, as we calculate full SO properties associated with 
each density peak, a single FOF group is naturally allowed to 
contain multiple SO host halos; thus the bridging or premerg- 
ing of FOF groups does not affect the final halo catalogs.. 

3.5. Calculating Halo Properties and Merger Trees 

Typically, several properties of interest are generated for 
halo catalogs. Regardless of the quality of particle assignment 
in the halo finder, careful attention to halo property calcula- 
tion is essential for consistent, unbiased results. 

3.5.1. Halo Positions and Velocities 



For positions, [Knebe et al.| ( |201 1} demonstrated that halo 
finders which calculated halo locations based on the maxi- 
mum density peak were more accurate than FOF-based halo 
finders which use the averaged location of all halo particles. 
The reason for this may be simply understood: as particle 
density rapidly drops in the outer reaches of a halo, the corre- 
sponding dispersion of particle positions climbs precipitously. 
Consequently, rather than increasing the statistical accuracy 
of the halo center calculation, including the particles at the 
halo boundary actually reduces it. The highest accuracy is in- 
stead achieved when the expected Poisson error (o-.v/a/A^) is 
minimized. As our halo finder has access (via the hierarchy 
of FOF subgroups) to the inner regions of the halo density dis- 
tribution, an accurate calculation of the center is possible by 
averaging the particle locations for the inner subgroup which 
best minimizes the Poisson error Typically, for a 10^ parti- 
cle halo, this estimator ends up averaging the positions of the 
innermost 10^ particles. 

The picture for halo velocities is not quite as simple. As 
noted in ^ the halo core can have a substantial velocity offset 
from the halo bulk. Since the galaxy hosted by the halo will 
presumably best track the halo core, we calculate the main 
velocity for the halo using the same Poisson-error minimiza- 
tion technique as for the halo positi on. F or calculating the 
bound/unbound mass of the halo (see |3.5.2"l ), however, we use 
the more appropriate averaged halo bulk velocity including 
substructure. 

3.5.2. Halo Masses 

For halo masses, Rockstar calculates spherical overdensi- 
ties according to a use r-specified density thresho ld: either the 
virial threshold, from Bryan & Norman ( 1998| l, or a density 
threshold relative to the background or to the critical density. 
As is usual, these overdensities are calculated using all the 
particles for all the substructure contained in a halo. On the 
other hand, subhalo masses have traditionally been a major 
point of ambiguity (for density-space halo finders). With a 
phase-space halo finder, such as Rockstar, the particles be- 
longing to the subhalo can be more reliably isolated from the 
host, and thus less ambiguity exists: the same method of cal- 
culating spherical overdensities may be applied to just the par- 
ticles belonging to the subhalo. In terms of the definition of 
where the subhalo "ends," Eq. l3]implies that the subhalo edge 
is effectively where the distribution of its particles in phase 
space becomes equidistant from the subhalo and its host halo. 



3.5.3. Unbinding Particles 

By default, Rockstar performs an unbinding procedure be- 
fore calculating halo mass and v„,fl^, although this may be 
switched off for studies of e.g., tidal remnants. Because the 
algorithm operates in phase space, the vast majority of halo 
particles assigned to central halos are actually b ound . We find 
typical boundedness values of 98% at z = 0; see 5 4.5 Even for 
substructure, unbound particles typically correspond to tidal 
streams at the outskirts of the subhalo, making a complicated 
unbinding algorithm unnecessary. For this reason, as well as 
to improve consistency of halo masses across timesteps, we 
perform only a single-pass unbinding procedure using a mod- 
ified Barnes-Hut method to accurately calculate particle po- 
tentials (see Appendix |b] for details)|j Since many users will 
be interested in classical halo finding only as opposed to re- 
covering tidal streams, the code by default does not output 
halos where fewer than 50% of the particles are bound; this 
threshold is user-adjustable, but changing it does not produce 
statistically significant effects on the clustering or mass func- 
tion until halo s wi th a bound fraction of less than 15% are 
included (see |4.5| l. 

We note that, in major mergers, a more careful approach to 
unbinding must be used. In many cases where merging halos 
initially have large velocity offsets, particles on the outskirts 
of the halos can mix in phase space before the halo cores 
themselves merge. This results in many particles being un- 
bound with respect to either of the two halo cores, even though 
they are bound to the overall merging system. As such, a naive 
unbinding of the particles would lead to the merging halos' 
masses dronping sharply for several timesteps prior to the fi- 
nal mergerPlTo counter this effect in major mergers, Rockstar 
calculates the gravitational potential of the combined merg- 
ing system, rather than for individual halos, when determining 
whether to unbind particles. This behavior is triggered by de- 
fault when two halos have a merger ratio of 1:3 or larger; this 
value is user-adjustable, but has little effect on the recovered 
mass function or clustering (see §4.5 i. 



3.5.4. Additional Halo Properties and Merger Trees 



the 



max^ 

max IS 
T. 



Two more common outputs of halo finders are v, 
maximum circular velocity and the scale radius, v, 
simply taken as the maximum of the quantity ■\/ GM(r)r 
it should be noted that this quantity is robust even for the 
smallest halos because of the extremely shallow dependence 
of v,„ax on radius. Our halo finder does not calculate R, ex- 
p hcitly from the halo p article distribution. Instead, we note as 
in Kly pin et al.| ( |20I1 ) that knowledge of both Vmax and M,,,-,- 
uniquely determmes the NFW ( [Navarro et al.|[T997] l profile. 
In particular, for NFW profiles, the radius iJ„,ai at which v,„ax 
is measured is a constant multiple b of the radius R,, and is 
given by: 

b-^ln(l + b) ^ 



d 
db 



1 + b 



= 



(5) 



This may be solved numerically, with the result that 

R,„,, = 2.1626R, (6) 

Instead of using R„,„^ directly (which is ill-determined for 
small halos), we make use of the ratio of Rmax/Rs to relate 

Provided enough interest, we may add the option of a multi-pass unbind- 
ing procedure in future versions of the halo finder. 

^ With the BDM halo finder jKlypin et al.|201l} , for example, we have 
observed halo masses which drop by a factor ot ttiree on account of this affect. 
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the mass enclosed within 2.1626/?j to v^^x, M^ir, and the con- 
centration C = Ryjr/Rs'- 



2.1626 



f(c) 



""GMyi, /(2.1626) 

/W = ln(l+x)--^ 
1 -\-x 



(7) 
(8) 



Thus, by numerically inverting the function on the left-hand 
side of Eq. |7j c may be found as a function of M,.,> and v,„„v 

We additionally calculate the angular momentum of the 
halo (using bound particles out to the desired halo radius) and 
the halo spin parameter (A), as introduced by Peebles ( 1969 ): 



A = 



GM: 



2.5 



(9) 



where J is the magnitude of the halo angular momentum and 
E is the total energy of the halo (potential and kinetic). 

Finally, we mention that our halo finder automatically cre- 
ates particle-based merger trees. For a given halo, its descen- 
dant is assigned as the halo in the next timestep which has the 
maximum number of particles in common (excluding parti- 
cles from subhalos). While it is possible to use these merger 
trees directly, we recommend inst ead to use the advanc ed 
merger tree algorithm discussed in |Behroozi et al.| ( 2011| l to 
further improve consistency across outputs of the halo finder 

4. TESTS & COMPARISONS 

The Rockstar algorithm has already underg one extensive 
testing and comparison to other halo finders in 'Kneb e et al.| 
(pOl 1 ). In tests with generated mock halos, Rockstar success- 
fully recovered halo properties for halos down to 20 particles, 
in many cases (e.g., for halo mass and bulk velocity) better 
than all seventeen other participating halo finders. In cases 
where it did not perform best, it was often only marginally 
behind one of the other phase-space halo finders. Notably, out 
of all the halo finders, it was the only one to fully successfully 
recover all halo properties (mass, location, position, velocity, 
and v,„ax) for a su bhalo coinciding wit h the center of its host 
halo. In addition, Knebe et al. (2011 1 compared mass func- 
tions, Vmax functions, correlation functions (for r > 2 Mpc), 
and halo bulk velocities for a cosmological simulation with 
1024^ particles; Rockstar had results comparable to other halo 
finders in all these results, although only the other phase-space 
halo finders shared its low mass completeness limit (^ 25 par- 
ticles forMaooc)- 

We t hus focus on compa risons beyond those already cov- 
ered in Knebe et al.| ( 201l] ). Our comparisons cover results 
for sever al di fferent dark matter simulations, briefly summa- 
rized in ^4.1 1 We provid e a v isual demonstration of the algo- 
rithm's performance in 5 4.2 a detailed comparison with the 
mass and correlation functions for other halo finders in ^4.3[ 
an ev aluation of the dynamical accuracy of halo properties in 
§4.4 and we pr esen t justification for our choice of the default 
parameters in ^4.5 1 Finally, we show figures demonstrating 
the excellent performance of the halo finder in ^4.6| 



4.1. Simulation Parameters 

We use five sets of simulations for this paper The first 
of these, Bolshoi, h as already been extensively detailed in 
|Klypin et al.| ( |20lT] l. To summarize the relevant parameters, 
it is a 2048'' particle simulation of comoving side length 250 
Mpc run using the ART simulation code (Kravtsov et al. 



1997| l on the NASA Ames Pleiades supercluster The as- 
sumed cosmology is f7„, = 0.27, Oa = 0.73, h = 0.7, = 0.96, 
and a = 0.82, similar to WMAP7 results ( |Komatsu et al.| 
201 1 1; the effective force-softening length is 1 kpc h~\ and 
the particle mass is 1.36 x lO^M© /i"'. 

We also use four simulations of different sizes from the 
LasDamas project (McBride et al, in preparation)]^ These 
have 1120^* to 1400-^ particles in comoving regions from 420 
Mpc /;"' to 2400 Mpc on a side, and were run using the 
GADGET-2 code (Springel 2005). Again, the assumed cos- 
mology is similar to WMAP7, with il,,, = 0.25, I^a = 0.75, 
h = 0.7, ris = 1.0, and a = 0.8. The effective force-softening 
lengths range from 8 to 53 kpc /i"', and the particle masses 
range from 1.87 x lO'^M© /i"' to 4.57 x lO^/i-'M©. Details 
of all the simulations are summarized in Table [T] 

4.2. Visual Demonstration 

In order to demonstrate how the algorithm performs on ha- 
los in a real simulation. Fig. [2]shows an example of how par- 
ticles are assigned to halos m a major merger; this example 
has been taken from the Bolshoi simulation at z = 0. From the 
top image in the figure, one might expect that two large ha- 
los separated by about 200 kpc are merging together, but 
careful analysis reveals that the larger halo actually consists 
of another major merger wherein the halo cores are separated 
by only 15 kpc /i"'. As shown in the bottom-left panel of 
the figure, the existence of the third massive halo is visible at 
a moderate significance level in position space — however, a 
position-space halo finder would have no way to correctly as- 
sign particles beyond the immediate locality of the two cores. 
Yet, in the bottom-right panel, the separation of the two halo 
cores is clearly distinguishable in velocity space. As such, not 
only can the distinct existence of the close-to-merging halos 
be reliably confirmed, but particle assignment to the two halos 
based on particle velocity coordinates can be reliably carried 
out as well. 

We also remark that phase-space halo-finding helps im- 
prove the accuracy of subhalo shapes by removing the need 
to perform a position-space cut to distinguish host particles 
from satellite particles. Satellites are usually offset in veloc- 
ity space from their hosts, but just as importantly, they usually 
also have a much lower velocity dispersion than their hosts. 
This implies that satellites may be reliably found even in the 
dense cores of halos — although the position space density is 
very high for host particles, the average velocity dispersion is 
large enough that the lower-dispersion subhalo particles can 
be reliably distinguished from host particles. Consequently, 
as shown in the middle-right panel of Fig. l2j satellites are 
picked out even in the dense halo centers witnout bias as re- 
gard to shape or size. 

4.3. Mass and Correlation Function Comparisons 

An extensive comparison of the mass function and corre- 
lation function between Rockstar and other h alo finders has 
already been conducted in |Knebe et al.| ( [20TT] ). Our algorithm 
has changed somewhat since that paper was published, and 
the mass and v„ma function comparison in that paper did not 
separate out the effects of central halos from satellites. In ad- 
dition, the correlation function comparison only compared the 
outskirts of halos beyond r = 2 Mpc and not smaller scales 
dominated by substructure. 
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TABLE 1 
Simulation Parameters 



Name 


Particles 


Size 


Particle Mass 


Softening 




Ua 


h 


o"8 




Code 


Bolshoi 


2048^ 


250 h- 


Mpc 


1.36 X lO'^h- 


Mq 


1 Ir' kpc 


0.27 


0.73 


0.7 


0.82 


0.96 


ART 


Consuelo 


1400^ 


420 r' 


Mpc 


1.87 X lO'/j- 


Me 


8 h-' kpc 


0.25 


0.75 


0.7 


0.8 


1.0 


GADGET-2 


Esmeralda 


1250^ 


640 r' 


Mpc 


9.31 X lO'/i- 




15 r' kpc 


0.25 


0.75 


0.7 


0.8 


1.0 


GADGET-2 


Carmen 


1120^ 


1000 h- 


' Mpc 


4.94 X lO'^r 




25 r' kpc 


0.25 


0.75 


0.7 


0.8 


1.0 


GADGET-2 


Oriana 


1280^ 


2400 h- 


' Mpc 


4.57 X 10"r 


'Me 


53 h-' kpc 


0.25 


0.75 


0.7 


0.8 


1.0 


GADGET-2 



Note. — Descriptions are in ^4.1| 

In Fig. [3] we present a comparison of the mass function 
for central halos in the Bolshoi and Las Damas simul ations. 
The re sults from Rockstar agree with those from Tink er et al.| 
(2008 1 at all masses for halos with more than 100 particles on 
the level of a few percent, well within the calibration speci- 
fication (5%) except when Poisson errors dominate. Detailed 
comparisons of the mass function from all of the Las Damas 
boxes, including fifty times the volume shown here, will be 
presented by McBride et al (in preparation). In the Bolshoi 
simulation, the mass function for centr al halos returned by 
Rockstar is virtually identical to BDM (Klypin et al. 201 l| l, 
differing by at most 5% over the entire mass range until Pois- 
son errors dominate; these differences stem mainly from dif- 
ferent conventions for unbinding particles, especially in major 
mergers. 

Figs. [4] and [5] show comparisons between the Vmax function 
for all halos and for satellites and the correlation function 
for Vmax-selected halos between R ockstar and the BD M al- 
gorithm for the Bolshoi simulation (jKlypin et al.| 201lll. The 
Vm-dx function for all halos is also mostly identical, differing 
by only 5% on average for halos above 100 km . For satel- 
lites, the deviations are slightly more significant, especially 
for very massive halos. Most notably, as a phase-space halo 
finder, Rockstar is able to track subhalos into the extreme in- 
ner reaches of halos. This enables it to track very massive 
subhalos even when their cores are very close to their hosts 
in position space, as long as they are sufficiently separated 
in velocity space. This is evident both in the increased num- 
ber density of subhalos with large v,nax as compared to BDM 
(Fig. PI, but also in the increased 1-halo contribution to the 
correlation function for massive halos (Fig.lSll. We also show 
comparisons for a small reg ion of Bolshoi ( 1/1 25th of the to- 
tal volume) where Subfind (^S pringel et al.|[2001i ) halos were 
also available in Fig. [5] For the halos considered in this latter 
comparison (Vmax > 100 km s"'), BDM may overpredict the 
number of major mergers within 20 kpc /?"', whereas Subfind 
may underpredict the number of major mergers within 30 kpc 
/i"', given the deviations seen in the correlation function as 
compared to larger scales. 

Fig.|6]shows the relationship between Mvir and Vmax for both 
satellites and centrals using the Rockstar halo finder on Bol- 
shoi; as may be expected, satellites have a very similar re- 
lation as compared to centrals. Due to dynamical stripping, 
however, the satellite relation has more scatter and is biased 
towards lower halo masses at fixed Vmax- 

4.4. Dynamical Tests 

In| Knebe et al.|p011| l, the authors performed a series of 



tests on mock halos m order to assess the accuracy of halo 
property recovery. While Rockstar performed extremely well 
in these tests, they nonetheless are not representative of the 
halos which one would expect to find in a real simulation (the 
tests considered only spherical, NFW/Plummer halos with 



TABLE 2 

Summary of algorithm parameters 



Variable Default Description 



Section 



h 0.28 Friends-of-friends linking length. 
/ 0.7 Fraction of particles in each hierarchical 
refinement level. 
A(z) virial Spherical overdensity definition 
BP 1 (Bound Properties) Whether halo prop- 
erties are calculated only using bound 
particles. 

UT 0.5 (Unbound Threshold) The minimum ra- 
tio of bound to total mass required for 
halos. 

MP 0.3 (Major-merger Potential) Mass ratio for 
mergers at which particles are evaluated 
for boundedness relative to the merg- 
ing system, as opposed to individual 
halos/subhalos. 



{3 



{73 



33^3 



T53 



333 



very little substructure). This is understandable — in a real 
simulation, there is no a priori "correct" answer for the re- 
covery of halo properties. Nonetheless, it is still possible to 
assess the precision of halo property recovery in a simulation 

by comparing halo properties between timest eps. 

To this end, we have used code from |Behroozi et al.|p011| l 
to check the halo property consistency for two position-space 



halo fin ders (BDM and Subfind: [Klypin et al.|201iySpruigeI 
|et al.|2001) and Rockstar on the Bolshoi simulation. This code 
simulates the gravitational evolution of halos purely based on 
their recovered properties (mass, scale radius, position, and 
velocity). Given the positions and velocities at one timestep, 
one can thus predict their positions and velocities at the next 
timestep with high accuracy; the difference between the pre- 
dicted and actual positions and velocities is one measure of 
the uncertainty in halo property recovery for the halo finder 

Fig. |7] shows the results of this analysis. Both BDM 
and Rockstar recover halo positions very precisely across 
timesteps; although Rockstar recovers positions better by a 
factor of ^ 2 for lower halo masses, both perform close to the 
force (softening) resolution of the simulation. Subfind, where 
halo positions are averaged over all particles, performs poorly 
in comparison. Rockstar recovers halo velocities much bet- 
ter (2-3 times as precisely) as compared to BDM, due largely 
to the fact that BDM uses only the inne rmost 100 particles 
to compute halo velocities (| Klypin et al.||201 1 ). By compari- 
son, Subfind appears to do extremely well; however, because 
it averages particle velocities over the entire halo instead of 
over the central region, it suffers from substantial systematic 
errors (as discussed in ^ which make its actual performance 
in recovering estimates of galaxy velocities always worse than 
Rockstar 

4.5. Evaluation of Default Parameters 
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Fig. 2. — Rockstar allows recovery of even very close major mergers. This figure shows an example of a major merger involving IO'^Mq halos from the 
Bolshoi simulation (Klypin et al. 2011). The top panel shows the complete particle distribution around the merging halos. In the second row, the left panel 
shows the host particle distribution, and the right panel shows the subhalo particle distribution, with particles colored according to subhalo membership. (The 
particle plotting size has been increased to show more clearly the extent of the small substructures in the right-hand panel). The two different colors in the 
left-hand panel hint at the fact that there are indeed three halos involved in the major merger, two of which are extremely close to merging. The uniform subhalo 
shapes in the right-hand panel suggest that subhalo particles can be distinguished without bias despite extreme variations in the host particle density between 
the subhalo centers and the subhalo outskirts. The bottom row shows more clearly the extremely close major merger. The bottom left-hand panel shows the 
full particle distribution in position space in a small region close to the merging halo cores; here, the bimodal distribution is evident, but distinguishing particle 
membership is impossible beyond the immediate vicinity of the cores. On the other hand, the bottom right-hand panel shows the same particles in velocity space, 
where the bimodality from the two cores shows a clear velocity separation, allowing particles to be reasonably assigned even in the halo bulk. 
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Halo Mass [M ] Halo Mass [M ] 



Fig. 3. — The halo mass function for distinct halos is very similar to previously published results. The top row shows comparisons to the Tinker et al. 1 2008J 
central halo mass function for each of the four LasDamas boxes (see Table[TJ. Good agreement is seen above 100 particles. The bottom row shows a comparison 
between the Rockstar and BDM (Klypin et al. 2011 1 halo finders on the Bolshoi simulation (2048^ particles, 250 Mpc /i"'). The left-han d plots show the full 
mass functions, and the right-hand plots show the residuals, with Poisson errors shown for the Bolshoi simulation. As noted in |Tinker et al.|j2008j , the calibrated 
mass range does n ot extend below 10' * br'MQ ; furthermore, the authors state "the behavior of the fitting function at lower masses is arbitrary." We therefore do 
not extrapolate the |Tinker et al.|j2008^ mass function in our comparisons. 




1000 



Fig. 4. — The halo velocity function is also very similar to previously published results. This figure shows comparisons between velocity (Vmax) functions for 
the Rockstar halo finder and BDM on the Bolshoi simulation (2048^ particles, 250 Mpc h^'). The left-hand plot shows all halos; the right-hand plot shows 
satellites only. 
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FlG. 5. — Two-point correlation functions are very similar to previously- 
published results except in the extreme centers of halos. This figure shows 
congelation functions for the Rockstar halo finder and BDM on the Bolshoi 
simulation (2048^ particles, 250 Mpc /j"'). The top panel shows the coiTela- 
tion function for Vmax > 300 km s"' (10,000 halos), the middle panel shows 
the correlation function for Vmax 

> 150 km (100,000 halos), and the bot- 
tom panel shows the correlation function for I'max > 100 km i"' . The bottom 
panel includes a comparison with Subfind. Subfind halos were only available 
for l/125th of Bolshoi (a total of ~ 1000 halos for this v^ax threshold). For 
a fair comparison, both Rockstar and BDM results for the bottom panel were 
computed on the same region of Bolshoi. 



Fig. 6. — Relationship between Myj, and Vn,ax for satellites (conditional 
density plot) and centrals (red line; error bars show the l-cr scatter about 
the mean) at z = in Bolshoi for the Rockstar halo finder While satellites 
have a very similar Mvir-i'max relation to centrals, the relation for satellites 
has slightly more scatter and is slightly biased towards lower masses at fixed 
Vmax, a consequence of dynamical stripping. 

In Fig. [HI we show the residual mass and correlation func- 
tions for changes in the default parameter settings for Rock- 
star All mass and correlation functions were calculated at 
z = for the Consuelo simulation. A summary of the main 
tunable parameters is shown in Table |2] 

In the top row of figures, it is evident that choosing a base 
3D linking length of b = 0.2 does not capture particles all the 
way out to the viiial radius for massive halos. Larger values of 
b (0.25-0.35) result in very similar mass functions; the stan- 
dard value of ^7 = 0.28 thus allows some degree of safety even 
for specific resimulations with unusual halo shapes. No effect 
on Vmax results from choosing larger linking lengths, and thus 
the correlation function is unaffected. For smaller values of b, 
FOF groups become more fragmented; at low particle num- 
bers, this can result in a single halo at high b being detected 
as multiple halos for a lower b. 

In the second row, the effects of varying the refinement 
threshold for the 6D FOF hierarchy are shown. The default 
behavior is to retain 70% of particles from the next higher re- 
finement level. If one retains more particles, then one is more 
likely to find low-significance objects — however, one is also 
more likely to find coincidental particle groups which do not 
correspond to actual halos. If one had a 90% particle retain- 
ment threshold, then approximately 4% of 30-particle halos 
(1O"M0 for Consuelo) would be false positives; this is ex- 
actly the residual difference between the mass functions for a 
90% threshold vs. the standard threshold in Fig.|8] However, a 
50% threshold is somewhat too low, as several percent fewer 
halos are detected as compared with the standard threshold. 
For halos with masses above 100 particles, where percent- 
level comparisons are more trustworthy (Tinker et al. 2008), 
the mass functions are almost indistinguishable. Some effect 
is seen in the autocoiTelation function for Vmax > 240 km s"' 
halos within 40 kpc. However, ^suggests that Rockstar with 
default parameters is already substantially complete for major 
mergers down to 10 kpc; as such, the additional halos found 
with / = 0.9 as opposed to the default / = 0.7 may represent 
false detections. 

In the third row, the effects of calculating halo properties 
using all particles (as opposed to the default behavior of only 
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Fig. 7. — Rockstar shows superior recovery of halo positions and velocities 
for cosmological halos as compared to both the BDM and Subfind halo find- 
ers. This figure shows a comparison of the joint consistency of halo position 
and velocity for two timesteps of the Bolshoi simulation at z = separated by 
40 Myr, including both distinct halos and subhalos. By comparing the evolu- 
tion of halo positions and velocities across timesteps to the values predicted 
by the laws of inertia and gravity, one may o btain an estimate of the reliabil- 
ity of halo property recovery (see text, also Behroozi at al.|201 I t. The Y-axis 
shows the difference between the predicted and measured position (top) and 
velocity (bottom). Rockstar offers excellent performance in locating halo po- 
sition centers, in all cases very close to the force resolution of the simulation; 
BDM performs almost as well, but Subfind peii'orms somewhat worse on ac- 
count of using a different position estimator. Rockstar offers significantly 
better performance than BDM for estimating halo velocities. Subfind's halo 
velocities appear more consistent than Rockstar, but because they represent 
the bulk velocity of the halo as opposed to the core velocity (which would cor- 
respond to measurable galaxy properties), they suffer from large systematic 
errors which make them significantly less accurate than Rockstar (see ^6). 
As such, the direct Subfind results are marked with an asterisk to signify 
caution in their interpretation. For reference, the systematic errors calculated 
in ^are included for Subfind in the green dash-double-dotted line ("Subfind 
w/ syst."). Subfind halos were only available for a small region of Bolshoi 
(l/125th of the total volume), and so error measurements do not extend to 
halo masses above lO'^'Mm. 



using bound particles) are shown. When using all particles, 
halo masses increase by 1-2% on average, implying that 98% 
of the initially-assigned particles were bound. Including un- 
bound particles affects satellite halos more than centrals, as 
evidenced by an increase in the correlation function in the 
regime where the 1-halo term is dominant. 

In the fourth row, the effects of changing the boundedness 
threshold requirement are shown. Lowering the threshold 
substantially (to UT=0.1) does not result in any more halos 
being found; however, raising the threshold to UT=0.7 results 
in fewer satellites being found and a sUght decrement in the 



mass function. 

In the fifth row, the effects of changing the major- merger 
potential threshold are shown; they have extremely little 
(<1%) effect on the mass function, and only a minor effect 
on the correlation function due to somewhat higher satellite 
masses for lower values of the threshold. 

4.6. Performance 

Figure [9] shows the excellent performance scaling of Rock- 
star The costliest part of our algorithm is the process of cal- 
culating the hierarchical refinement levels, which takes more 
time for halos with more substructure. For that reason, struc- 
ture finding takes more time both for higher resolution simu- 
lations and for late redshifts (the runtime per snapshot scales 
approximately as T cx cP '^^). Nonetheless, the halo finder is 
still so efficient (5-10 CPU hours per billion particles at z = 0) 
that high-throughput access to particle data is important to 
avoid wasting CPU time. Indeed, despite its advanced capa- 
bilities, it is roughly 4-5 times as efficient as other halo find- 
ers that include substructure with published performance data 
(jKnollmann & Knebe 2009; Springel et al. 2010), as shown 
in Fig. [9] Some caution is necessary in comparing these tim- 
ings directly, as the systems used were not the same; however, 
the processors used for calculating Rockstar's performance 
are over four years old, negating any advantage from newer 
technol ogy. The most mod em friends-of-friends halo finders 
(e.g., Wood ring et al.|2()Tl ) are faster by roughly a factor of 
two than Rockstar ori large datasets; however, these cannot 
identify substructure nor can they produce complete catalogs 
of halos with SO masses. 

5. ESTIMATING SATELLITE COMPLETENESS LIMITS 

For many decades, early computational simulations suf- 
fered from the so-called "overmerging" problem, where satel- 
lite halos were found to disappear almost as soon as they 
passed wit hin the boundary of a larger cluster (see Klypin] 
et al.|1999 for a discussion). This situation has improved dra- 
matically on account of increased mass and force resolution in 
simulations; however, as recently as for the Millennium Sim- 
ulation, simulators have found that they needed to add "or- 
phan" satellite galaxies (galaxies which exist without a corre- 
sponding subhalo) in order to match the small-scale clustering 
seen in observations ( |Kitzbichler & W hite 2008 ). With recent 
observations pushing the previous small-sc ale limit of cluster- 
ing measurement s from 100 kpc to 10 kpc ( [Tinker et al. 201 1 
Jiang et al.|2011| , it is important to understand the limitations 
of cosmological simulations' abilities to recover substructure 
at this level. 

In observations, the distributions of satellite galaxies 
around such massive hosts is completely self-similar (density 
proportional to a pow er law of radius, wi th an exponent be- 
tween -1.7 and -1.5, [Tink er et al.||2011|l. However, Figure 
[T0| acutely shows the tension between this observational re- 
sult and what is found in the Bolshoi simulation for cluster- 
sized halos (7?vh ^ 1 Mpc); in the Bolshoi simulation, there 
is a clear radial incompleteness scale in the recovered volume 
density for merger ratios below 1:30. 

Making a quantitative estimate of the radial incompleteness 
scale is difficult because the true satellite distribution, /3(r), is 
unknown. We can nonetheles s connect the luminosity con- 
straints in I Tinker et al.| pOl l| l to constraints on dark matter 
halos via the assumption that galaxy luminosity is tightly cor- 
related with halo peak mass (i.e., the largest mass that pro- 
genitors of a given halo or subhalo have ever had). This as- 
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Fig. 8. — Dependence of the mass function and correlation functions (Vm„v > 460 km s"' and i',„a,- > 240 km s"') for halos in Consuelo (1400^ particles, 
420 Mpc /i"') on the parameter choices in the Rockstar halo finder (Table|2l. Mass function error bars are Poisson, and a 1O"M0 halo has approximately 30 
particles; coiTelation function eiTor bars are jackknife. The top row showsnow the mass function recovery is affected by changing the linking length used to 
preselect bound groups; clearly, h = 0.2 is too small, but the default value b = 0.28 gives results indistinguishable from larger values. The second row shows 
how the mass function is affected when changing the fraction of particles retained in each successive 6D refinement level. The default fraction (0.7) represents 
the optimal balance between recovery of small halos while at the same time being resilient to noise at low particle number or high particle density (see text). 
The third row shows that roughly 1-2% of halo particles are unbound; as shown in the correlation function plot on the right-hand side, satellites tend to have 
inore unbound particles than centrals. The fourth row shows that raising the threshold for the required fraction of bound particles mainly affects completeness 
for close halo pairs; otherwise, changing the threshold has no effect. The fifth row shows alinost negligible effects on the mass function and coiTelation function 
from including additional contributions to the gravitational potential; the main effect is better seen in merger trees. 
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Fig. 9. — Rockstar demonstrates extremely efficient use of CPU time as 
compared to other halo finders with published results. Figure shows how 
CPU time requirements per billion particles analyzed at z = depend on sim- 
ulation particle mass for a number of simulations (Bolshoi and LasDamas 
simulations, as well as a lightcone for the DBS Blind Cosmology Challenge). 
Rockstar tests were timed on four-year-old 2.3 GHz AMD Opter on 2376 
cores. Comparisons are include d to two other halo finders ( Knoll mann &| 
[Knebe 2009 Spring el et al.|2010 > with published performance data (note that 
mese were not run on the same system, and so may not be directly compara- 
ble). For Rockstar, typical runtime depends on scale factor approximately as 
T oc a^'^; at earlier times, fewer particles are clustered and analysis is more 
efficient. 

sumption has been tested in a large number of studies and 
been found to accurately reproduce both the halo mass to stel- 
lar mass relation and the luminosity-dependent clustering of 
galaxies CBehroozi et al."2010', 'Moster et al.'2010[|Guo et al. 



2010; Conroy & Wechsler 2009; Wang & Jing 2010[ |Conroy 
et al.| |2006). As concerns halos, this assumption would im- 



ply that the radial dependence of satellites with a given peak 
mass should also follow a power law with an exponent be- 
tween -1.7 and -1.5. 

We thus conservatively define the satellite incompleteness 
scale as the radius at which the logarithmic slope of the satel- 
lite density p{r) becomes shallower than -1.5; steeper num- 
bers could overestimate the true logarithmic slope of the satel- 
lite radial distribution. However, we include the possibility 
that the true cutoff slope should be -1.7 in our treatment of 
the errors. We show the results for the radial completeness 
scale as a function of satellite mass ratio (satellite peak mass 
compared to host mass) for massive hosts (lO'^-^M,^ < M < 
10''*Mq) in Fig. 
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We also include results for the BDM 
halo finder in Fig.|10|to show the comparison with non-phase- 
space halo finding^ne benefit of phase-space finding is strik- 
ing for satellite mass ratios above 1:30, with Rockstar able to 
find major mergers consistently down to a tiny fraction of the 
virial radius. However, for satellites with smaller mass ratios, 
Rockstar is unable to perform any better than BDM. This find- 



ing is especially unexpected considering the results in Knebe 
|et al.| ( [20lT| l, wherein Rockstar was the only halo finder able 
to accurately recover all the halo properties of a satellite (mass 
ratio 1:100) placed directly at the center of a million-particle 
IO'^Mq halo. 

In order to explain such a drastic reduction in ability to find 
satellites, it is instructive to consider the effects of tidal strip- 
ping on a satellite in a simulation with limited force and mass 
resolution. In particular, due to gravitational softening, the 
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Fig. 10. — Top panel: A plot showing the radial distribution of existing 
satellite halos at z = in cluster-scale halos as a function of the satellite mass 
ratio (p, the ratio between the peak satellite mass and the host mass) for halos 
found with Rockstar in the Bolshoi simulation. Bottom panel: the radial 
completeness limit as a function of satellite mass ratio for cluster-scale halos 
in Bolshoi for both Rockstar and the BDM halo finder, defined as the radius 
at which the logarithmic slope of the satellite distribution reaches -1.5 (see 
text). Error bars include Poisson uncertainties as well as unce rtainties in thi s 
cutoff slope, which could reasonably be as high as -1.7 (Tink er et al.|201 1) . 
The upturn in radial incompleteness for major mergers in BDM is a conscious 
choice of BDM's author. A theoretical limit is also shown, but the calculation 
is valid only for the assumption of circular orbits and spherical NFW profiles, 
and thus carries large systematic errors. The asterisks thus denote that some 
caution is necessary in interpreting these results. (See discussion in 01. 



maximum density of a halo cannot continue increasing indef- 
initely towards its center; instead, it will threshold at some 
distance fies from the center For halos with a low enough 
mass, there may not be any particles within f^es of the center; 
for these halos, the effective maximum density is degraded 
even further Due to this effect, satellites are vulnerable to 
tidal disruption much earlier for larger values of /„.5 and for 
larger particle masses. In particular, the fluid Roche limit (un- 
der the assumption of a significantly more massive host) is 
given by 



t^Roche ~ 2.4/?H 



(10) 



where liRoche is the radius within which tidal disruption occurs, 
Rh is the host radius, ph is the enclosed averaged host density. 
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and ps is the average satellite density. In terms of halos, this 
means that tidal disruption will occur at a distance R from the 
host center where the enclosed host halo average density is 
2.4"^ ~ 7% of the satellite peak density ps- 

A naive application of the Roche limit under the assumption 
of a spherical NFW profile, circular satellite orbits, the mean 
mass-concentration relation from Bolshoi, as well as the mass 
and force resolutions of Bolshoi yields an estimate of the the- 
oretical disintegration limit shown in Figure 10 A number 



of significant limitations apply to this calculation, because it 
does not include any mass stripping for satellites after they en- 
ter the host, and it does not include any effects from eccentric 
orbits, where satellites at a given radius would be expected 
to have passed through closer radii many times. The latter 
bias is stronger for low-mass satellites, which experience less 
dynamical friction and so are expected to orbit many more 
times before destruction than larger halos; as such, low-mass 
satellites at a given radius are on average more stripped and 
have had more exposure to high tidal fields than high-mass 
satellites at the same radius. This calculation also does not 
include numerical simulation artifacts wherein forces on one 
part of a satellite may be calculated differently from forces 
on another part (e.g., direct summation versus multipole ex- 
pansion). Nonetheless, many of these biases go in the direc- 
tion of earlier destruction of satellites; as such, the calculation 
represents a useful lower limit on the radius where satellites 
may exist in massive clusters. For a much more detailed as- 



sessm ent of subhalo completeness, we refer the reader to Wu 
|etaL| ( |20lT| i. 

This calculation and the results in Fig. [TO] would suggest 
that, regardless of the halo finding technique used, simulations 
of clusters still suffer from an "overmerging" problem at the 
very innermost radii due to limited force and mass resolution. 
Depending on the scientific questions being answered with a 
simulation, these missing/incomplete halos may be more or 
less relevant. For a high resolution simulation like Bolshoi, 
the incompleteness does not impact the correlation function 
on scales larger than ~' 100 kpc, and is present only in the 
inner regions of massive halos. In general, when determin- 
ing the desired mass and force resolution to required, the re- 
quirements to resolve satellites in the inner regions of massive 
clusters will be much more stringent. In these regions, "or- 
phan" galaxies may represent an attractive solution to reduce 
computational requirements, but the need for these galaxies is 
somewhat mitigated with increased force resolution and with 
the effective halo finding in inner regions available with Rock- 
star. 

6. HALO CORE VELOCITY OFFSETS 

In the ACDM paradigm, the main method by which halos 
reach a relaxed state is dynamical friction; namely, energy in 
bulk motions relative to the halo center of gravity is trans- 
formed into increased halo velocity dispersion. Dynamical 
friction depends not only on the background density that a 
satellite halo is passing through, but also on the satellite mass 
and velocity. For a massive host halo, the high background 
density means that incoming satellites will initially transfer 
momentum to the host with high efficiency. However, mas- 
sive host halos also have very strong tidal fields, and once 
satellites are disrupted into cold velocity streams, dynamical 
friction becomes much less effective at transferring momen- 
tum into the inner halo regions. This combination of dynam- 
ical friction and tidal forces suggests that momentum transfer 
may be more efficient in the outer regions of a host halo and 
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Fig. 1 1 . — Variation between halo core velocity and halo bulk velocity as 
a function of halo mass and redshift. Panels show median absolute offsets 
between the halo bulk velocity and radially-binned average particle velocity 
as a function of redshift and halo mass. Calculations including substructure 
("All Part.") as well as excluding particles in substructure ("Excl. SS") ai'e 
shown. All panels show central halos from the Bolshoi simulation; the mid- 
dle and lower panels also show central halos from the Consuelo simulation. 
In all cases, each radial bin contains at least 400 particles so as to minimize 
the effects of Poisson noise. For the lower panel, the numbers of halos used 
in Bolshoi were 18, 226, and 484 at redshifts z = 1.41, 0.53, and 0.00, respec- 
tively. All other panels (as well as Consuelo for z = 0.00) averaged results 
from more than 2000 halos. Each panel shows halos in a mass bin of 0.48 
dex (1 to 3 times IO'^Mq, IO'^^Mq, and IO'^'Mq, respectively); the respec- 
tive virial radii are ~ 0.3 /i-'Mpc, ~ 0.6 /r'Mpc, and ~ 1.3 /j-'Mpc. 
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Fig. 12. — Velocity offsets of the particle velocities at a given radius from 
the halo bulk velocity. There is no distance at which they match consistently. 
Nonetheless, inner radial bins (excluding substructure) often have consistent 
velocities, suggesting that the halo core velocity is both well-defined and 
physical. The top two panels show the .v-velocity offsets between the halo 
bulk velocity and radially-binned halo particle velocity for 10 central halos 
with masses between 10'** and 3 X IO'^Mq randomly selected from the Bol- 
shoi simulation at z = 0. The bottom panel shows z-velocity offsets for the 
same halos. The top panel shows calculations excluding substructure; the 
middle and bottom panels include all particles. Matching colors between 
the panels correspond to the same halo. Radial bins are evenly spaced in 
log(r), and as such contain significantly more particles at large radii. The 
red- violet halo (dashed line) is undergoing a major merger in its outskirts, as 
such, the halo bulk velocity is over 400 km s"' offset along the A"-direction 
from the halo core velocity. 



Fig. 13. — There is a large intrinsic scatter in the core-bulk velocity off- 
set, at least ±0.3 dex among halos at a given mass. These panels show the 
conditional probability density of offsets between the halo bulk velocity and 
radially-binned halo particle velocity for 484 central halos with masses be- 
tween lO'** and 3 X IQ''*Mq from the Bolshoi simulation at z = 0. The red 
line shows the median of the offsets, and the error bars encompass it Icr of 
the distribution. 
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Poisson error minimization method) and the bulk velocity of all halo particles 
within Svir as a function of halo mass at z = 0. This is also compared to the 
halo particle velocity dispersion, ct, . 
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less efficient in the inner regions, leading to an offset between 
the mean velocity of the central density peak and the bulk ve- 
locity of the halo. 

We examine the presence of this effect by calculating 
radially-averaged halo particle velocities (i.e., particle veloci- 
ties averaged in spherical shells) and comparing to the halo 
bulk velocity (averaged over all halo particles) for a large 
number of central halos in both the Bolshoi and Consuelo 
simulations over a range of halo masses (10'^ to IO^^Mq) and 
redshifts (z = to 1.4). We additionally consider the effects 
of excluding particles belonging to substructure from calcula- 
tions of the radially-binned velocities; results for all of the 
median offsets (both including and excluding substructure) 



are shown in Fig. 11 

Figure [TT] demonstrates that the core-bulk velocity differ- 
ence can bequite dramatic at high redshift for massive halos: 
halos with 10"*Mq < M < 3 x IO'^^Mq at z = 1.4 have a me- 
dian offset of over 450 km s"' between the velocity averaged 
within 0.1/?vir (excluding substructure) and the bulk velocity 
of the halo. This effect is lower for both smaller halos and 
later times, presumably due to lower velocity offsets in in- 
coming substructure for smaller halos and a reduced merger 
rate at later times. Nonetheless, a clear difference between 
the inner and bulk velocities is present in all halos tested even 
down to z = 0. In all cases, there is a definite radial trend, with 
the velocity in the inner region of the halo being farther away 
from the bulk velocity than the halo outskirts, consistent with 
the hypothesis of reduced efficiency of momentum transfer in 
the inner regions of the halo. 

In the calculations which exclude substructure, a clear 
plateau is evident below 0.1-0.27?vir where the core veloc- 
ity offset stabilizes, at least for halos with M > IQ^^Mq; we 
have insufficient mass resolution to test this effect robustly 
for smaller halos. This transition is extremely consistent 
with expectations for where most of the satellite mass will 
be stripped in such halos. For a satellite halo with a con- 
centration of c = 10 falling into a massive host with c = 5, 
Eq. 10 would suggest that 90% of the satellite mass (and 



momentum) will be stripped from the main satellite density 
peak at radii greater than O.37?vir,host, with another 9% stripped 
by 0.17?vir.ho,st- As tidal forces increase dramatically towards 
the center of a halo, they should disrupt efficient momentum 
transfer within 0.1/?vir.host even for satellites on highly radial 
orbits. 

In the calculations which include substructure, the median 
velocity offsets are less than the results which exclude sub- 
structure because substructure is included when calculating 
the bulk halo velocity. In the very innermost regions of ha- 
los, however, both radial velocity averages converge, suggest- 
ing that our results should be robust to the subhalo finding 
technique used. On the outskirts of halos, the presence of 
non-disrupted subhalos results in an upturn in the offset be- 
tween the radially-averaged velocity (including substructure) 
and the bulk halo velocity. 

Due to significantly higher mass resolution, the Bolshoi 
simulation is able to probe the velocity offsets to significantly 
smaller radii as compared to the Consuelo simulation. Fur- 
thermore, the Bolshoi simulation is able to better resolve sub- 
structure than Consuelo; as such, the difference between ra- 
dial velocity offset profiles which include and exclude sub- 
structure in Consuelo is less than the difference for Bolshoi. 
Except for this caveat, the median velocity offsets are in ex- 
cellent agreement between the two simulations. 



One more interesting feature of Fig. 1 1 is that the median 
velocity offset never reaches zero for any of the halos in either 
simulation. That is to say, there is no single radius where the 
bulk-averaged halo velocity corresponds to the actual average 
motion of particles in a radial shell. We examine this issue fur- 
ther by considering ten halos randomly drawn from the Bol- 
shoi simulation at z = in the mass range lO'"* to 3 x IQ^'^Mq. 
As before, we compute the average velocity in radial bins for 
each halo, but instead of the absolute offset from the bulk ve- 
locity, we show just the offsets between the x-components of 
the two velocities in Fig. [12] 

The results for individualhalos in Fig. [T2]are just as striking 
as for the median offsets in Fig. 1 1 Almost every halo in Fig. 



11 exhibits the same velocity offset plateau within O.IR^ 
suggesting that the velocity within 0.1/? vii has a physical in- 
terpretation corresponding to the actual bulk motion of those 
particles. If substructure is excluded, the radially-averaged x- 
velocity often never matches the bulk jc- velocity, implying an 
overall offset between the main halo motion and satellites. If 
substructure is included, the radially-averaged x-velocity can 
briefly match the bulk jc-velocity (and it must do so, by the 
mean value theorem), but there is no reason why the radially 
averaged velocities along other axes must match the bulk ve- 
locity at the same time; in fact, as shown for the z-velocities 
in the bottom panel of Fig. [121 they do not. As such, the bulk 
velocity of the particles in thehalos shown here does not actu- 
ally correspond to the velocity averaged over any of the radial 
shells. 

To investigate this latter point in more detail, we show con- 
ditional density plots of the absolute offsets between radially- 
averaged velocities and the bulk velocity for the complete 
sample of 10'"* to 3 x 10'"* Mq (again at z = in Bolshoi) in 
Figure 13 While there is significant scatter in the absolute 
offsets (on the order of 0.3 dex, after correcting for the vari- 
ation over the bin width), only a tiny fraction (<C 1 %) of the 
halos ever have a radial bin in which the averaged velocity 
matches the bulk velocity to better than 20 km s"', and most 
have significant offsets (^ 100 km s"') at all radii from the 
halo center. 

These res ults su pport the use of the Poisson minimization 
method of 5 3.5.1 as a velocity estimator for Rockstar, as the 
method primarily probes the inner regions of halos with ve- 

Indeed, for both 



12 



locity profiles like those shown in Fig. 
the Bolshoi and Consuelo simulations, we find an extremely 
consistent me an vel ocity offset between the core velocity (as 
discussed in § 3.5.1[ ) and the halo bulk velocity, as shown in 
Fig. 14 Significantly, these offsets are much larger than either 
t he dire ct Poisson error estimates in determining velocities 
( P.5.1| l or the cross-timestep velocity error estimates in Fig. 
[7|( ^4.4| i. The main differences between Bolshoi and Consuelo 
come from reduced mass resolution and substructure resolu- 
tion in the latter simulation, resulting in a velocity estimate 
slightly closer to the bulk average for Consuelo as compared 
to Bolshoi. 

These velocity offsets display a power-law-like behavior 
across more than four orders of magnitude in halo mass; how- 
ever, this power law is not the same as the power law for 
the average velocity dispersion ct,, (which is proportional to 
M'/^); instead, the velocity offset relative to the velocity dis- 
persion increases for increasing halo masses. At z = 0, the 
average core-velocity offset may be fit by the empirical rela- 
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tion: 



(Ay2) =4kms- 



M 



0.46 



(11) 



This relation implies that, for the most massive M ^ 2 x 
1O'^M0 clusters at z = 0, the average core-bulk velocity offset 
is of order 400 km s"' . 



7. CONCLUSIONS 

We have presented a novel method of finding halos in phase 
space, the Rockstar halo finder. This method has several key 
strengths which, in combination, improve on previous ap- 
proaches to robustly identifying halos and their substructures: 

1. high acc uracy of recovered halo properties (as dis- 
cussed in |Knebe et al.|20l"T] l, 

2. the abi lity to consistently and usefully define a subhalo 
mass ( p^^ , 

3. the ability to reliably track major mergers down to the 
inner host halo core (©land 



4. explicit grid-independence, orientation-independence, 
and halo shape-independence for recovered halos (Q, 

5. calibrated estimates of the accuracy of recovered halo 
properties in full cosmological simulations (^4.4i, 



6. extremely efficient resource consumption both in CPU 
time and memory (§4.6 1, 



7. massively parallel implementation, supporting large fu- 
ture simulations (Appendix [A|, and 

8. a publicly available codebase: 

http: //code . google .com/p/rockstar I ■ 

In addition, the ha lo finder complements the method of 
|Behroozi et al.| ( [20TT i) for creating merger trees, making it part 
of the first halo finding system to recover halos using seven 
dimensions of information (phase space plus time). These 
properties make it ideal for a large number of scientific appli- 
cations, especially for studies of subhalos and subhalo particle 
distributions. We have demonstrated the ability to probe sub- 
structure into the very inner regions of halos, at least as far as 
simulations are capable of resolving. This is critical for com- 
parison with galaxy statistics in dense regimes, and paves the 
way for direct comparisons of clustering between simulations 
and observations in the regions where the unknown effects of 
baryons are the most significant. 

We have also investigated the radial dependence of halo ve- 
locity. We find that there exists a well-defined core veloc- 
ity (for r < O.lRyir) which is always distinct from the bulk- 
averaged halo velocity. Using the bulk velocity to estimate 
galaxy peculiar velocities thus carries significant systematic 
error, which can be up to 400 km s"' in the most massive 
clusters at z = and is larger at higher redshifts. Furthermore 
(as discussed in the halo velocity averaged in radial bins 
never matches thebulk average velocity; as such, the bulk av- 
erage velocity is only an ensemble property of the entire halo, 
and the local dynamics have a much more complicated and in- 
teresting structure than can be described by the bulk velocity 
and the velocity dispersion alone. 
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APPENDIX 
LOAD BALANCING 

The simplest approach to load balancing — equal volumes — works well on scales where the universe is homogenous, namely, 
when individual analysis regions are at least 100 Mpc on a side. However, cosmic variance on smaller scales requires a more 
sophisticated approach for algorithms designed to run on small simulations or very large numbers of processors. 

The approach we use is to divide the simulation volume into chunks requiring approximately equal memory. We divide the 
number of available processors, A'^, into three factors (A^:c> Ny, and A'^), which control the number of divisions along each of the 
principal axes of the simulation. Naturally, A^i, N^, and A^j. are chosen to be as close as possible to N^^^ so as to minimize the 
surface-to-volume ratio of the divisions. Initially, particles are binned into B = 10* bins according to their x-coordinate, and the 
simulation is divided into A^^^ regions of equal particle count along the x-axis. Then, within each of these regions, particles are 
binned according to their >'-coordinate, and each region is split accordingly along the 3'-axis into Ny subregions of equal particle 
count. Finally, each subregion is split in the same way into sub- subregions along the z-axis. This process is accurate and 
efficient, requiring a fixed total data transfer size of 0{BN^I^), an amount which is completely independent of the total number 
of particles. 

3D FOF group calculation is performed in parallel on these regions without any inter-process communication. However, for 
phase-space analysis, the number of recovered groups can vary considerably depending on the overall density of the analysis 
region (e.g., voids will have lower fractions of bound particles). For that reason, FOFs are distributed to processors in small work 
units; as soon as a processor finishes one work unit, it receives another to work on (possibly from a completely different part 
of the simulation), so that maximum concurrency can be maintained at all times. However, a single FOF group is never (in the 
current implementation) split up for analysis among multiple processors. Thus, the analysis time for the largest FOF group sets 
the lower hmit for the wall clock time taken to analyze the simulation snapshot, regardless of the number of processors available 
for use. 



20 



BEHROOZI, WECHSLER, & WU 



UNBINDING 



To calculate which particles are bound to halos, we use a modified Barnes-Hut algorithm (|Barn es & jjutl |1986| l to calculate 
particle potential energies with a binary space partitioning (BSP) tree|^ As noted in Salmon & Warren ( 199417 the cell refinement 



criterion suggested by Barnes and Hut gives unacceptably high maximum errors (as opposed to average errors); in addition, our 
use of a BSP tree (which gives fewer wasted refinement levels than an octree) renders the Barnes and Hut criterion somewhat 
inapplicable. Instead, given a cell containing n particles with mass center xq, we calculate a threshold monopole approximation 
distance — that is, the distance from xq at which a monopole wo uld be an appropria te approximation of the potential. We can 



estimate the relative error 9 in the calculated potential following Salmon & Warren] ( | 199 4) as a function of distance d from the 
center xq as being bounded by: 

^W<Tir^ (Bi) 

^.^EkHk^^ (B2) 

where x, and m,- are the locations and masses of the n particles in the cell, and r,nax is the maximum distance between xq and 
any of the particles in the cell. Thus, for a given choice of the relative error 0, this equation may be inverted to give a minimum 
distance for acceptable use of the monopole approximation: 



^(e)>^ + /("^j (B3) 

This equation is valid regardless of the boundary sizes of the cell; hence it is appropriate for use even with a BSP tree. We choose 
9 to give potentials accurate to 4%; from our tests, this is sufficient to correctly estimate the boundedness of 99.8% of all halo 
particles. 



^max ' ^ ^ ^ 



' Several other halo finders, such as BDM and AHF, use the assumption 
of spherical symmetry to calculate particle potentials; we find that, while this 



works well for the vast majority of halos, the potential calculated in this way 
can be dramatically wrong for halos undergoing major mergers, especially if 
the halo center is incorrectly identified. 



